function [Tea, varargout] = alphaVAR(alphat,Q,K,q,V_prior,T_prior,stat,sv)
% Calculate (V)AR coefficients of state equation
% input:
% alphat is TxK and Q KxK

% output: T is K x K*p

% Prepare input for multiReg
Y = alphat(q+1:end,:);

X =[];
for qx=q:-1:1
    X = [X,alphat(qx:end-q+qx-1,:)];
end

if sv>1
    Sigma=Q;
else
    if sv==1
        Sigma = exp(Q);
    end
    
    if sv==0
        Sigma = exp(Q).^2;
    end
    
end

Teaf = multiReg(Y,X,Sigma,V_prior,T_prior,sv);

TeaMult = reshape(Teaf,[K*q,K])';

Phi_statCheck = VARcompanion(TeaMult,eye(K),0);

% Stationarity check
inx =1;
stato=0;
if stat==1
    while (max(abs( eig(Phi_statCheck))) >= 1);
        Teaf = multiReg(Y,X,Sigma,V_prior,T_prior,sv);
        
        TeaMult = reshape(Teaf,[K*q,K])';
        Phi_statCheck = VARcompanion(TeaMult,eye(K),0);
        inx=inx+1;
        
        if inx>10000
            disp('Problems with stationarity');
            stato =1;            
          %  break;          
        end
    end
end
%% Set zero coefficients from prior to zero
Tea = zeros(K,K*q);


Tea(1:K,1:K) = diag(diag(TeaMult(1:K,1:K)));

for qx=2:q
    Tea(1:K,1+K*(qx-1):K*(qx-1)+K) = diag(diag(TeaMult(1:K,1+K*(qx-1):K*(qx-1)+K)));
end

if sv==1 % Watanabe and Omori stochastic volatility sampler
    nu_star = Y - X*Tea';
else % for Kim, Shephard, and Chib
    nu_star = Y - X*Tea';
    nu_star=log(nu_star.^2+0.001)';
end

varargout(1) = {stato};
varargout(2) = {nu_star};


